home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_200 / 247_01 / bnround.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-04-17  |  4.0 KB  |  188 lines
  1. /*
  2.  *   MIRACL euclidean mediant rounding routine
  3.  *   bnround.c
  4.  */
  5.  
  6. #include <stdio.h>
  7. #include "miracl.h"
  8. #define abs(x)  ((x)<0? (-(x)) : (x))
  9.  
  10. /* Access global variables */
  11.  
  12. extern int  depth;    /* error tracing .. */
  13. extern int  trace[];  /* .. mechanism     */
  14. extern int  nib;      /* length of bigs   */
  15. extern small base;    /* number base      */
  16. extern bool check;    /* overflow check   */
  17.  
  18. extern big w0;        /* workspace variables */
  19. extern big w5;
  20. extern big w6;
  21. extern big w7;
  22.  
  23. int euclid(x,num)
  24. big x;
  25. int num;
  26. { /* outputs next c.f. quotient from gcd(w5,w6) */
  27.     static long u,v;
  28.     static small a,b,c,d;
  29.     static int   n,oldn,q;
  30.     static bool last,carryon;
  31.     small r,m;
  32.     long lr,lq;
  33.     big t;
  34.     if (num==0)
  35.     {
  36.         oldn=(-1);
  37.         carryon=FALSE;
  38.         last=FALSE;
  39.         if (compare(w6,w5)>0)
  40.         { /* ensure w5>w6 */
  41.             t=w5,w5=w6,w6=t;
  42.             return (q=0);
  43.         }
  44.     }
  45.     else if (num==oldn || q<0) return q;
  46.     oldn=num;
  47.     if (carryon) goto middle;
  48. start:
  49.     if (size(w6)==0) return (q=(-1));
  50.     n=w5[0];
  51.     carryon=TRUE;
  52.     a=1;
  53.     b=0;
  54.     c=0;
  55.     d=1;
  56.     if (n==1)
  57.     {
  58.         last=TRUE;
  59.         u=w5[1];
  60.         v=w6[1];
  61.     }
  62.     else
  63.     {
  64.         m=w5[n]+1;
  65.         if (sizeof(long)==2*sizeof(small))
  66.         { /* use longs if they are double length */
  67.             if (n>2)
  68.             { /* squeeze out as much significance as possible */
  69.                 u=muldiv(w5[n],base,w5[n-1],m,&r);
  70.                 u=u*base+muldiv(r,base,w5[n-2],m,&r);
  71.                 v=muldiv(w6[n],base,w6[n-1],m,&r);
  72.                 v=v*base+muldiv(r,base,w6[n-2],m,&r);
  73.             }
  74.             else
  75.             {
  76.                 u=(long)base*w5[n]+w5[n-1];
  77.                 v=(long)base*w6[n]+w6[n-1];
  78.                 last=TRUE;
  79.             }
  80.         }
  81.         else
  82.         {
  83.             u=muldiv(w5[n],base,w5[n-1],m,&r);
  84.             v=muldiv(w6[n],base,w6[n-1],m,&r);
  85.         }
  86.     }
  87. middle:
  88.     forever
  89.     { /* work only with most significant piece */
  90.         if (last)
  91.         {
  92.             if (v==0) return (q=(-1));
  93.             lq=u/v;
  94.         }
  95.         else
  96.         {
  97.             if (((v+c)==0) || ((v+d)==0)) break;
  98.             lq=(u+a)/(v+c);
  99.             if (lq!=(u+b)/(v+d)) break;
  100.         }
  101.         if (lq>=TOOBIG || lq>=MAXBASE/abs(d)) break;
  102.         q=lq;
  103.         r=a-q*c;
  104.         a=c;
  105.         c=r;
  106.         r=b-q*d;
  107.         b=d;
  108.         d=r;
  109.         lr=u-lq*v;
  110.         u=v;
  111.         v=lr;
  112.         return q;
  113.     }
  114.     carryon=FALSE;
  115.     if (b==0)
  116.     { /* update w5 and w6 */
  117.         check=OFF;
  118.         divide(w5,w6,w7);
  119.         check=ON;
  120.         if (lent(w7)>nib) return (q=(-2));
  121.         t=w5,w5=w6,w6=t;   /* swap(w5,w6) */
  122.         copy(w7,x);
  123.         return (q=size(x));
  124.     }
  125.     else
  126.     {
  127.         check=OFF;
  128.         premult(w5,c,w7);
  129.         premult(w5,a,w5);
  130.         premult(w6,b,w0);
  131.         premult(w6,d,w6);
  132.         add(w5,w0,w5);
  133.         add(w6,w7,w6);
  134.         check=ON;
  135.     }
  136.     goto start;
  137. }
  138.  
  139.  
  140. void round(num,den,z)
  141. big num,den;
  142. flash z;
  143. { /* reduces and rounds the fraction num/den into z */
  144.     int s;
  145.     if (ERNUM) return;
  146.     if (size(num)==0) 
  147.     {
  148.         zero(z);
  149.         return;
  150.     }
  151.     depth++;
  152.     trace[depth]=34;
  153.     if (TRACER) track();
  154.     if (size(den)==0)
  155.     {
  156.         berror(13);
  157.         depth--;
  158.         return;
  159.     }
  160.     copy(num,w5);
  161.     copy(den,w6);
  162.     s=exsign(w5)*exsign(w6);
  163.     insign(PLUS,w5);
  164.     insign(PLUS,w6);
  165.     if (compare(w5,w6)==0)
  166.     {
  167.         convert((small)s,z);
  168.         depth--;
  169.         return;
  170.     }
  171.     if (size(w6)==1)
  172.     {
  173.         if (w5[0]>nib)
  174.         {
  175.             berror(13);
  176.             depth--;
  177.             return;
  178.         }
  179.         copy(w5,z);
  180.         insign(s,z);
  181.         depth--;
  182.         return;
  183.     }
  184.     build(z,euclid);
  185.     insign(s,z);
  186.     depth--;
  187. }
  188.